CAMB Python example notebook

Run it online yourself in Binder.

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import sys, platform, os
import matplotlib
from matplotlib import pyplot as plt
import numpy as np
import camb
from camb import model, initialpower
print('Using CAMB %s installed at %s'%(camb.__version__,os.path.dirname(camb.__file__)))
# make sure the version and path is what you expect

Using CAMB 1.4.1 installed at c:\work\dist\git\camb\camb

#Set up a new set of parameters for CAMB
pars = camb.CAMBparams()
#This function sets up with one massive neutrino and helium set using BBN consistency
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.06, omk=0, tau=0.06)
pars.InitPower.set_params(As=2e-9, ns=0.965, r=0)
pars.set_for_lmax(2500, lens_potential_accuracy=0);

#calculate results for these parameters
results = camb.get_results(pars)

#get dictionary of CAMB power spectra
powers =results.get_cmb_power_spectra(pars, CMB_unit='muK')
for name in powers: print(name)

total
unlensed_scalar
unlensed_total
lensed_scalar
tensor
lens_potential

#plot the total lensed CMB power spectra versus unlensed, and fractional difference
totCL=powers['total']
unlensedCL=powers['unlensed_scalar']
print(totCL.shape)
#Python CL arrays are all zero based (starting at L=0), Note L=0,1 entries will be zero by default.
#The different CL are always in the order TT, EE, BB, TE (with BB=0 for unlensed scalar results).
ls = np.arange(totCL.shape[0])
fig, ax = plt.subplots(2,2, figsize = (12,12))
ax[0,0].plot(ls,totCL[:,0], color='k')
ax[0,0].plot(ls,unlensedCL[:,0], color='C2')
ax[0,0].set_title('TT')
ax[0,1].plot(ls[2:], 1-unlensedCL[2:,0]/totCL[2:,0]);
ax[0,1].set_title(r'Fractional TT lensing')
ax[1,0].plot(ls,totCL[:,1], color='k')
ax[1,0].plot(ls,unlensedCL[:,1], color='C2')
ax[1,0].set_title(r'$EE$')
ax[1,1].plot(ls,totCL[:,3], color='k')
ax[1,1].plot(ls,unlensedCL[:,3], color='C2')
ax[1,1].set_title(r'$TE$');
for ax in ax.reshape(-1): 
    ax.set_xlim([2,2500])
    ax.set_xlabel(r'$\ell$');

(2551, 4)

# The lensing B modes are non-linear, so need to be calculated carefully if you want them accurate (even at low ell)
# Need both high lmax, non-linear lensing and high k 
# lens_potential_accuracy=1 turns on the latter, and can be increased to check precision

pars.set_for_lmax(2500, lens_potential_accuracy=1)
results = camb.get_results(pars)
lmax2500CL = results.get_lensed_scalar_cls(CMB_unit='muK')
                  
pars.set_for_lmax(4000, lens_potential_accuracy=1)
results = camb.get_results(pars)
lmax4000CL = results.get_lensed_scalar_cls(CMB_unit='muK')

pars.set_for_lmax(4000, lens_potential_accuracy=2)
results = camb.get_results(pars)
accCL = results.get_lensed_scalar_cls(CMB_unit='muK')

pars.set_for_lmax(6000, lens_potential_accuracy=4)
results = camb.get_results(pars)
refCL = results.get_lensed_scalar_cls(CMB_unit='muK')

fig, ax = plt.subplots(1,2, figsize = (12,4))
ax[0].plot(ls,totCL[:len(ls),2], color='C0')
ax[0].plot(ls,lmax2500CL[:len(ls),2], color='C1')
ax[0].plot(ls,lmax4000CL[:len(ls),2], color='C2')
ax[0].plot(ls,accCL[:len(ls),2], color='C3')
ax[0].plot(ls,refCL[:len(ls),2], color='k')

ax[0].set_xlim([2,2500])
ax[0].set_xlabel(r'$\ell$',fontsize=13)
ax[0].set_ylabel(r'$\ell(\ell+1)C_\ell^{BB}/2\pi\,[\mu {\rm K}^2]$', fontsize=13)

ax[1].plot(ls[2:],totCL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C0')
ax[1].plot(ls[2:],lmax2500CL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C1')
ax[1].plot(ls[2:],lmax4000CL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C2')
ax[1].plot(ls[2:],accCL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C3')

ax[1].axhline(0,color='k')
ax[1].set_xlim([2,2500])
ax[1].set_xlabel(r'$\ell$',fontsize=13)
ax[1].set_ylabel('fractional error', fontsize=13);
ax[1].legend(['Default accuracy','lmax=2500, lens_potential_accuracy=1',
           'lmax=4000, lens_potential_accuracy=1','lmax=4000, lens_potential_accuracy=2']);

#You can calculate spectra for different primordial power spectra without recalculating everything
#for example, let's plot the BB spectra as a function of r
pars.set_for_lmax(4000, lens_potential_accuracy=1)
pars.WantTensors = True
results = camb.get_transfer_functions(pars)
lmax=2000
rs = np.linspace(0,0.2,6)
for r in rs:
    inflation_params = initialpower.InitialPowerLaw()
    inflation_params.set_params(ns=0.96, r=r)
    results.power_spectra_from_transfer(inflation_params) #warning OK here, not changing scalars
    cl = results.get_total_cls(lmax, CMB_unit='muK')
    plt.loglog(np.arange(lmax+1),cl[:,2])
plt.xlim([2,lmax])
plt.legend(["$r = %s$"%r for r in  rs], loc='lower right');
plt.ylabel(r'$\ell(\ell+1)C_\ell^{BB}/ (2\pi \mu{\rm K}^2)$')
plt.xlabel(r'$\ell$');

WARNING:root:power_spectra_from_transfer with non-linear lensing does not recalculate the non-linear correction

#Now get matter power spectra and sigma8 at redshift 0 and 0.8
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(ns=0.965)
#Note non-linear corrections couples to smaller scales than you want
pars.set_matter_power(redshifts=[0., 0.8], kmax=2.0)

#Linear spectra
pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)
s8 = np.array(results.get_sigma8())

#Non-Linear spectra (Halofit)
pars.NonLinear = model.NonLinear_both
results.calc_power_spectra(pars)
kh_nonlin, z_nonlin, pk_nonlin = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)

Note: redshifts have been re-sorted (earliest first)

print(results.get_sigma8())

[0.53304915 0.80260802]

for i, (redshift, line) in enumerate(zip(z,['-','--'])):
    plt.loglog(kh, pk[i,:], color='k', ls = line)
    plt.loglog(kh_nonlin, pk_nonlin[i,:], color='r', ls = line)
plt.xlabel('k/h Mpc');
plt.legend(['linear','non-linear'], loc='lower left');
plt.title('Matter power at z=%s and z= %s'%tuple(z));

#Plot CMB lensing potential power for various values of w
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(As=2e-9, ns=0.965)
pars.set_for_lmax(2000, lens_potential_accuracy=1)

ws = np.linspace(-1.5, -0.6, 5)
for w in ws:
    pars.set_dark_energy(w=w, wa=0, dark_energy_model='fluid') 
    results = camb.get_results(pars)
    cl = results.get_lens_potential_cls(lmax=2000)
    plt.loglog(np.arange(2001), cl[:,0])

plt.legend(ws)
plt.ylabel('$[L(L+1)]^2C_L^{\phi\phi}/2\pi$')
plt.xlabel('$L$')
plt.xlim([2,2000]);

---

You can view the parameters (as used by fortran CAMB internals) by just printing the parameter object. See the docs for parameter and structure descriptions

# parameters can also be read from text .ini files, for example this sets up a best-fit 
# Planck 2018 LCDM cosmology (base_plikHM_TTTEEE_lowl_lowE_lensing). 
# [Use planck_2018_acc.ini if you need high-ell and/or accurate BB and CMB lensng spectra at beyond-Planck accuracy]
pars = camb.read_ini('https://raw.githubusercontent.com/cmbant/CAMB/master/inifiles/planck_2018.ini')
# for a local github installation you can just do 
# pars=camb.read_ini(os.path.join(camb_path,'inifiles','planck_2018.ini'))
print(pars)

class: <CAMBparams>
 WantCls = True
 WantTransfer = True
 WantScalars = True
 WantTensors = False
 WantVectors = False
 WantDerivedParameters = True
 Want_cl_2D_array = True
 Want_CMB = True
 Want_CMB_lensing = True
 DoLensing = True
 NonLinear = NonLinear_both
 Transfer: <TransferParams>
   high_precision = True
   accurate_massive_neutrinos = False
   kmax = 1.3464234
   k_per_logint = 0
   PK_num_redshifts = 1
   PK_redshifts = [0.0]
 want_zstar = False
 want_zdrag = False
 min_l = 2
 max_l = 2700
 max_l_tensor = 600
 max_eta_k = 18000.0
 max_eta_k_tensor = 1200.0
 ombh2 = 0.0223828
 omch2 = 0.1201075
 omk = 0.0
 omnuh2 = 0.0006451439
 H0 = 67.32117
 TCMB = 2.7255
 YHe = 0.2454006
 num_nu_massless = 2.046
 num_nu_massive = 1
 nu_mass_eigenstates = 1
 share_delta_neff = True
 nu_mass_degeneracies = [0.0]
 nu_mass_fractions = [1.0]
 nu_mass_numbers = [1]
 InitPower: <InitialPowerLaw>
   tensor_parameterization = tensor_param_indeptilt
   ns = 0.9660499
   nrun = 0.0
   nrunrun = 0.0
   nt = 0.0
   ntrun = 0.0
   r = 0.0
   pivot_scalar = 0.05
   pivot_tensor = 0.05
   As = 2.100549e-09
   At = 0.0
 Recomb: <Recfast>
   min_a_evolve_Tm = 0.0011098779505118728
   RECFAST_fudge = 1.125
   RECFAST_fudge_He = 0.86
   RECFAST_Heswitch = 6
   RECFAST_Hswitch = True
   AGauss1 = -0.14
   AGauss2 = 0.079
   zGauss1 = 7.28
   zGauss2 = 6.73
   wGauss1 = 0.18
   wGauss2 = 0.33
 Reion: <TanhReionization>
   Reionization = True
   use_optical_depth = True
   redshift = 10.0
   optical_depth = 0.05430842
   delta_redshift = 0.5
   fraction = -1.0
   include_helium_fullreion = True
   helium_redshift = 3.5
   helium_delta_redshift = 0.5
   helium_redshiftstart = 6.0
   tau_solve_accuracy_boost = 1.0
   timestep_boost = 1.0
   max_redshift = 50.0
 DarkEnergy: <DarkEnergyPPF>
   w = -1.0
   wa = 0.0
   cs2 = 1.0
   use_tabulated_w = False
 NonLinearModel: <Halofit>
   Min_kh_nonlinear = 0.005
   halofit_version = mead2016
   HMCode_A_baryon = 3.13
   HMCode_eta_baryon = 0.603
   HMCode_logT_AGN = 7.8
 Accuracy: <AccuracyParams>
   AccuracyBoost = 1.0
   lSampleBoost = 1.0
   lAccuracyBoost = 1.0
   AccuratePolarization = True
   AccurateBB = False
   AccurateReionization = True
   TimeStepBoost = 1.0
   BackgroundTimeStepBoost = 1.0
   IntTolBoost = 1.0
   SourcekAccuracyBoost = 1.0
   IntkAccuracyBoost = 1.0
   TransferkBoost = 1.0
   NonFlatIntAccuracyBoost = 1.0
   BessIntBoost = 1.0
   LensingBoost = 1.0
   NonlinSourceBoost = 1.0
   BesselBoost = 1.0
   LimberBoost = 1.0
   SourceLimberBoost = 1.0
   KmaxBoost = 1.0
   neutrino_q_boost = 1.0
 SourceTerms: <SourceTermParams>
   limber_windows = True
   limber_phi_lmin = 100
   counts_density = True
   counts_redshift = True
   counts_lensing = False
   counts_velocity = True
   counts_radial = False
   counts_timedelay = True
   counts_ISW = True
   counts_potential = True
   counts_evolve = False
   line_phot_dipole = False
   line_phot_quadrupole = False
   line_basic = True
   line_distortions = True
   line_extra = False
   line_reionization = False
   use_21cm_mK = True
 z_outputs = []
 scalar_initial_condition = initial_adiabatic
 InitialConditionVector = []
 OutputNormalization = 1
 Alens = 1.0
 MassiveNuMethod = Nu_trunc
 DoLateRadTruncation = True
 Evolve_baryon_cs = False
 Evolve_delta_xe = False
 Evolve_delta_Ts = False
 Do21cm = False
 transfer_21cm_cl = False
 Log_lvalues = False
 use_cl_spline_template = True
 SourceWindows = []
 CustomSources: <CustomSources>
   num_custom_sources = 0
   c_source_func = None
   custom_source_ell_scales = []
 

#The dark energy model can be changed as in the previous example, or by assigning to pars.DarkEnergy.
#e.g. use the PPF model
from camb.dark_energy import DarkEnergyPPF, DarkEnergyFluid
pars.DarkEnergy = DarkEnergyPPF(w=-1.2, wa=0.2)
print('w, wa model parameters:\n\n', pars.DarkEnergy)
results = camb.get_background(pars)

#or can also use a w(a) numerical function 
#(note this will slow things down; make your own dark energy class in fortran for best performance)
a = np.logspace(-5, 0, 1000)
w = -1.2 + 0.2 * (1 - a)
pars.DarkEnergy= DarkEnergyPPF()
pars.DarkEnergy.set_w_a_table(a, w)
print('Table-interpolated parameters (w and wa are set to estimated values at 0):\n\n' 
      ,pars.DarkEnergy)
results2 = camb.get_background(pars)

rho, _ = results.get_dark_energy_rho_w(a)
rho2, _ = results2.get_dark_energy_rho_w(a)
plt.plot(a, rho, color='k')
plt.plot(a, rho2, color='r', ls='--')
plt.ylabel(r'$\rho/\rho_0$')
plt.xlabel('$a$')
plt.xlim(0,1)
plt.title('Dark enery density');

w, wa model parameters:

 class: <DarkEnergyPPF>
 w = -1.2
 wa = 0.2
 cs2 = 1.0
 use_tabulated_w = False
 
Table-interpolated parameters (w and wa are set to estimated values at 0):

 class: <DarkEnergyPPF>
 w = -1.2
 wa = 0.1993346429183863
 cs2 = 1.0
 use_tabulated_w = True
 

#Get various background functions and derived parameters
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
results = camb.get_background(pars)
print('Derived parameter dictionary: %s'%results.get_derived_params())

Derived parameter dictionary: {'age': 13.740451258557629, 'zstar': 1090.5875300814807, 'rstar': 144.2080769525876, 'thetastar': 1.0446246836476256, 'DAstar': 13.804774021712861, 'zdrag': 1059.2314028925762, 'rdrag': 146.98719294834376, 'kd': 0.14067153336071683, 'thetad': 0.1617760955140736, 'zeq': 3441.8515054149993, 'keq': 0.010503792217322724, 'thetaeq': 0.8079811517064835, 'thetarseq': 0.44697513328564176}

z = np.linspace(0,4,100)
DA = results.angular_diameter_distance(z)
plt.plot(z, DA)
plt.xlabel('$z$')
plt.ylabel(r'$D_A /\rm{Mpc}$')
plt.title('Angular diameter distance')
plt.ylim([0,2000])
plt.xlim([0,4]);

print('CosmoMC theta_MC parameter: %s'%results.cosmomc_theta())

CosmoMC theta_MC parameter: 0.010443924812710163

#You can also directly access some lower level quantities, for example the CMB transfer functions:
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
data = camb.get_transfer_functions(pars)
transfer = data.get_cmb_transfer_data()
print('Number of sources (T, E, phi..): %s; number of ell: %s; number of k: %s '%tuple(transfer.delta_p_l_k.shape))

Number of sources (T, E, phi..): 3; number of ell: 87; number of k: 2737 

#Plot the transfer functions as a function of k for various ell
fig, axs = plt.subplots(2,2, figsize=(12,8), sharex = True)
for ix, ax in zip([3, 20, 40, 60],axs.reshape(-1)):
    ax.plot(transfer.q,transfer.delta_p_l_k[0,ix,:])
    ax.set_title(r'$\ell = %s$'%transfer.L[ix])
    if ix>1: ax.set_xlabel(r'$k \rm{Mpc}$')

#Note that internal samplings can be quite sparsely sampled, e.g. look at l=2 transfer function
def plot_cmb_transfer_l(trans, ix):
    _, axs = plt.subplots(1,2, figsize=(12,6))
    for source_ix, (name, ax) in enumerate(zip(['T', 'E'], axs)):
        ax.semilogx(trans.q,trans.delta_p_l_k[source_ix,ix,:])
        ax.set_xlim([1e-5, 0.05])
        ax.set_xlabel(r'$k \rm{Mpc}$')
        ax.set_title(r'%s transfer function for $\ell = %s$'%(name, trans.L[ix]))
plot_cmb_transfer_l(transfer,0)

#So if you want to make nice plots, either interpolate or do things at higher than default accuracy
pars.set_accuracy(AccuracyBoost=2) #higher accuracy, so higher sampling density
data = camb.get_transfer_functions(pars)
plot_cmb_transfer_l(data.get_cmb_transfer_data(),0)
pars.set_accuracy(AccuracyBoost=1); #re-set default

#Similarly for tensor transfer function 
#e.g. see where various C_L^BB come from in k by plotting normalized transfer**2 (C_l is ~ integral over log k P(k) T(k)^2)
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.WantScalars = False
pars.WantTensors = True
pars.set_accuracy(AccuracyBoost=2)
data = camb.get_transfer_functions(pars)
transfer = data.get_cmb_transfer_data('tensor')
print(r'Calculated L: %s'%transfer.L)
plt.figure(figsize=(14,3))
ixs=[13,19,21]
ls = [transfer.L[i] for i in ixs]
cols=['b','r','c']
for ix,col in zip(ixs, cols):
    k_weight = transfer.delta_p_l_k[2,ix,:]**2
    k_weight /= np.sum(k_weight)
    plt.semilogx(transfer.q,k_weight, color=col)
plt.xlim([1e-3, 0.1])
plt.legend(ls)
plt.xlabel(r'$k \rm{Mpc}$')
plt.title(r'Contribution to B from primordial tensor power spectrum for various $\ell$')
#compare with k_* = l/chi*, note DAstar is in GPc, so multiply by 1000 to get standard Mpc units used for k
derived = data.get_derived_params()
for l,col in zip(ls,cols):
    plt.axvline(l/(1000*derived['DAstar']), color=col, ls=':', lw=2)

Calculated L: [  2   3   4   5   6   7   8   9  10  11  12  13  14  15  17  19  21  23
  25  27  29  31  33  35  37  40  45  50  55  60  65  70  75  80  85  90
 110 130 150 175 200 250 300 350 400 450 500 550 600]

#if you want to combine the transfer functions with the primordial power spectra, you can get the latter via
k=10**np.linspace(-5, 1, 50)
pars.InitPower.set_params(ns=0.96, r=0.2) #this functions imposes inflation consistency relation by default
scalar_pk= pars.scalar_power(k)
tensor_pk= pars.tensor_power(k)
plt.semilogx(k,scalar_pk);
plt.semilogx(k,tensor_pk);
plt.xlabel(r'$k \rm{Mpc}$')
plt.ylabel(r'${\cal P}(k)$')
plt.legend(['scalar', 'tensor']);

#set_params is a shortcut routine for setting many things at once
pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, As=2e-9, ns=0.95)
data= camb.get_background(pars)

#There are functions get plot evolution of variables, e.g. for the background as a function of conformal time:
eta = 10**(np.linspace(1, 4,300))
back_ev = data.get_background_time_evolution(eta, ['x_e', 'visibility'])
fig, axs= plt.subplots(1,2, figsize=(12,5))
axs[0].semilogx(eta, back_ev['x_e'])
axs[1].loglog(eta, back_ev['visibility'])
axs[0].set_xlabel(r'$\eta/\rm{Mpc}$')
axs[0].set_ylabel('$x_e$')
axs[1].set_xlabel(r'$\eta/\rm{Mpc}$')
axs[1].set_ylabel('Visibility');
fig.suptitle('Ionization history, including both hydrogen and helium recombination and reionization');

#or as a function of redshift
z = 10**np.linspace(2, 4, 300)
back_ev = data.get_background_redshift_evolution(z, ['x_e', 'visibility'], format='array')
fig, axs= plt.subplots(1,2, figsize=(12,5))
for i, (ax, label), in enumerate(zip(axs, ['$x_e$','Visibility'])):
    ax.semilogx(z, back_ev[:,i])
    ax.set_xlabel('$z$')
    ax.set_ylabel(label)
    ax.set_xlim([500,1e4])

# ..and perturbation transfer functions, e.g. for k=0.1. Note that quantities are synchronous gauge unless otherwise specified
# Also note v_newtonian_cdm, v_newtonian_baryon include a factor of -k/H, and Weyl a factor of k^2 -
# see https://camb.readthedocs.io/en/latest/transfer_variables.html
print('Available variables are %s'%camb.model.evolve_names)

Available variables are ['k/h', 'delta_cdm', 'delta_baryon', 'delta_photon', 'delta_neutrino', 'delta_nu', 'delta_tot', 'delta_nonu', 'delta_tot_de', 'Weyl', 'v_newtonian_cdm', 'v_newtonian_baryon', 'v_baryon_cdm', 'a', 'etak', 'H', 'growth', 'v_photon', 'pi_photon', 'E_2', 'v_neutrino', 'T_source', 'E_source', 'lens_potential_source']

eta = np.linspace(1, 400, 300)
ks = [0.02,0.1]
ev = data.get_time_evolution(ks, eta, ['delta_baryon','delta_photon'])
_, axs= plt.subplots(1,2, figsize=(12,5))
for i, ax in enumerate(axs):
    ax.plot(eta,ev[i,:, 0])
    ax.plot(eta,ev[i,:, 1])
    ax.set_title('$k= %s$'%ks[i])
    ax.set_xlabel(r'$\eta/\rm{Mpc}$');
plt.legend([r'$\Delta_b$', r'$\Delta_\gamma$'], loc = 'upper left');

#or as a function of redshift
z = np.linspace(500,5000,300)
ks = [0.02,0.1]
ev = data.get_redshift_evolution(ks, z, ['delta_baryon','delta_cdm', 'delta_photon'])
_, axs= plt.subplots(1,2, figsize=(12,5))
for i, ax in enumerate(axs):
    ax.plot(z,ev[i,:, 0])
    ax.plot(z,ev[i,:, 1])
    ax.plot(z,ev[i,:, 2])
    ax.set_title(r'$k= %s/\rm{Mpc}$'%ks[i])
    ax.set_xlabel('$z$');
plt.legend([r'$\Delta_b$', r'$\Delta_c$', r'$\Delta_\gamma$'], loc = 'upper right');

#Here you can see oscillation of delta_photon, subsequent decay of the potential and change to Mezsaroz growth in delta_cdm
eta = 10**(np.linspace(0, 3, 500))
def plot_ev(ev, k):
    plt.figure(figsize=(8,6))
    plt.loglog(eta,ev[:,0])
    plt.loglog(eta,np.abs(ev[:,1]))
    plt.loglog(eta,-ev[:,2]/k**2)  # Weyl is k^2*phi
    plt.title(r'$k= %s/\rm{Mpc}$'%k)
    plt.xlabel(r'$\eta/\rm{Mpc}$');
    plt.legend([r'$\Delta_c$', r'$|\Delta_\gamma|$', r'$-(\Phi+\Psi)/2$'], loc = 'upper left');

k=0.3
plot_ev(data.get_time_evolution(k, eta, ['delta_cdm','delta_photon', 'Weyl']),k)

#Note that time evolution can be visually quite sensitive to accuracy. By default it is boosted, but you can change this. e.g.
plot_ev(data.get_time_evolution(k, eta, ['delta_cdm','delta_photon', 'Weyl'],lAccuracyBoost=1),k)
plot_ev(data.get_time_evolution(k, eta, ['delta_cdm','delta_photon', 'Weyl'],lAccuracyBoost=10),k)

#It's also possible to plot quantities in other gauges, or arbitrary symbolic expressions,
#using camb.symbolic.
#For example, this plots the Newtonian gauge density compared to the synchronous gauge one
import camb.symbolic as cs
Delta_c_N = cs.make_frame_invariant(cs.Delta_c, 'Newtonian')
ev=data.get_time_evolution(k, eta, ['delta_cdm',Delta_c_N])
plt.figure(figsize=(6,4))
plt.loglog(eta,ev[:,0])
plt.loglog(eta,ev[:,1])
plt.title(r'$k= %s/\rm{Mpc}$'%k)
plt.xlabel(r'$\eta/\rm{Mpc}$');
plt.legend([r'$\Delta_c^{\rm synchronous}$', r'$\Delta_c^{\rm Newtonian}$'], fontsize=14);

# Or see that the Newtonian-gauge CDM peculiar velocity decays roughly propto 1/a on sub-horizon 
# scales during radiation domination after the potentials have decayed so there is no driving force
# (leading to logarithmic Meszaros growth of the CDM density perturbation)
k=4
vc_N = cs.make_frame_invariant(cs.v_c, 'Newtonian')
ev=data.get_time_evolution(k, eta, [Delta_c_N, vc_N, 'a'])
plt.figure(figsize=(6,4))
plt.loglog(eta,ev[:,0])
plt.loglog(eta,-ev[:,1])
eta_ln=6*np.sqrt(3)*np.pi/4/k
horizon_index =  np.searchsorted(eta, eta_ln, side='right')
plt.loglog(eta,-ev[horizon_index,1]*ev[horizon_index,2]/ev[:,2], ls='--')
plt.ylim([1e-4,5e2])
plt.title(r'$k= %s/\rm{Mpc}$'%k)
plt.xlabel(r'$\eta/\rm{Mpc}$');
plt.legend([r'$\Delta_c^{\rm Newtonian}$',r'$-v_c^{\rm Newtonian}$',r'$v_c^{\rm Newtonian}\propto 1/a$'], fontsize=14);

For further details of camb.symbolic and examples of what can be done see the CAMB symbolic ScalEqs notebook (run now in Binder). This also serves as documentation for the scalar equations implemented in CAMB.

#For calculating large-scale structure and lensing results yourself, get a power spectrum
#interpolation object. In this example we calculate the CMB lensing potential power
#spectrum using the Limber approximation, using PK=camb.get_matter_power_interpolator() function.
#calling PK(z, k) will then get power spectrum at any k and redshift z in range.

nz = 100 #number of steps to use for the radial/redshift integration
kmax=10  #kmax to use
#First set up parameters as usual
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(ns=0.965)

#For Limber result, want integration over \chi (comoving radial distance), from 0 to chi_*.
#so get background results to find chistar, set up a range in chi, and calculate corresponding redshifts
results= camb.get_background(pars)
chistar = results.conformal_time(0)- results.tau_maxvis
chis = np.linspace(0,chistar,nz)
zs=results.redshift_at_comoving_radial_distance(chis)
#Calculate array of delta_chi, and drop first and last points where things go singular
dchis = (chis[2:]-chis[:-2])/2
chis = chis[1:-1]
zs = zs[1:-1]

#Get the matter power spectrum interpolation object (based on RectBivariateSpline). 
#Here for lensing we want the power spectrum of the Weyl potential.
PK = camb.get_matter_power_interpolator(pars, nonlinear=True, 
    hubble_units=False, k_hunit=False, kmax=kmax,
    var1=model.Transfer_Weyl,var2=model.Transfer_Weyl, zmax=zs[-1])

#Have a look at interpolated power spectrum results for a range of redshifts
#Expect linear potentials to decay a bit when Lambda becomes important, and change from non-linear growth
plt.figure(figsize=(8,5))
k=np.exp(np.log(10)*np.linspace(-4,2,200))
zplot = [0, 0.5, 1, 4 ,20]
for z in zplot:
    plt.loglog(k, PK.P(z,k))
plt.xlim([1e-4,kmax])
plt.xlabel('k Mpc')
plt.ylabel('$P_\Psi\, Mpc^{-3}$')
plt.legend(['z=%s'%z for z in zplot]);

Now do integral to get convergence power spectrum, using Limber approximation ($k\approx (l+0.5)/\chi$) $$ C_l^\kappa \approx [l(l+1)]^2\int_0^{\chi_*} d\chi \left( \frac{\chi_*-\chi}{\chi^2\chi_*}\right)^2 P_\Psi\left(\frac{l+0.5}{\chi}, z(\chi)\right) $$ where $P_\Psi $ is obtained from the interpolator.

#Get lensing window function (flat universe)
win = ((chistar-chis)/(chis**2*chistar))**2
#Do integral over chi
ls = np.arange(2,2500+1, dtype=np.float64)
cl_kappa=np.zeros(ls.shape)
w = np.ones(chis.shape) #this is just used to set to zero k values out of range of interpolation
for i, l in enumerate(ls):
    k=(l+0.5)/chis
    w[:]=1
    w[k<1e-4]=0
    w[k>=kmax]=0
    cl_kappa[i] = np.dot(dchis, w*PK.P(zs, k, grid=False)*win/k**4)
cl_kappa*= (ls*(ls+1))**2

#Compare with CAMB's calculation:
#note that to get CAMB's internal calculation accurate at the 1% level at L~2000, 
#need lens_potential_accuracy=2. Increase to 4 for accurate match to the Limber calculation here
pars.set_for_lmax(2500,lens_potential_accuracy=2)
results = camb.get_results(pars)
cl_camb=results.get_lens_potential_cls(2500) 
#cl_camb[:,0] is phi x phi power spectrum (other columns are phi x T and phi x E)

#Make plot. Expect difference at very low-L from inaccuracy in Limber approximation, and
#very high L from differences in kmax (lens_potential_accuracy is only 2, though good by eye here)
cl_limber= 4*cl_kappa/2/np.pi #convert kappa power to [l(l+1)]^2C_phi/2pi (what cl_camb is)
plt.loglog(ls,cl_limber, color='b')
plt.loglog(np.arange(2,cl_camb[:,0].size),cl_camb[2:,0], color='r')
plt.xlim([1,2000])
plt.legend(['Limber','CAMB hybrid'])
plt.ylabel('$[L(L+1)]^2C_L^{\phi}/2\pi$')
plt.xlabel('$L$');

#The non-linear model can be changed like this:
pars.set_matter_power(redshifts=[0.], kmax=2.0)
pars.NonLinearModel.set_params(halofit_version='takahashi')
kh_nonlin, _, pk_takahashi = results.get_nonlinear_matter_power_spectrum(params=pars)
pars.NonLinearModel.set_params(halofit_version='mead')
kh_nonlin, _, pk_mead = results.get_nonlinear_matter_power_spectrum(params=pars)

fig, axs=plt.subplots(2,1, sharex=True, figsize=(8,8))
axs[0].loglog(kh_nonlin, pk_takahashi[0])
axs[0].loglog(kh_nonlin, pk_mead[0])
axs[1].semilogx(kh_nonlin, pk_mead[0]/pk_takahashi[0]-1)
axs[1].set_xlabel(r'$k/h\, \rm{Mpc}$')    
axs[1].legend(['Mead/Takahashi-1'], loc='upper left');

#Get angular power spectrum for galaxy number counts and lensing
from camb.sources import GaussianSourceWindow, SplinedSourceWindow

pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(As=2e-9, ns=0.965)
pars.set_for_lmax(lmax, lens_potential_accuracy=1)
#set Want_CMB to true if you also want CMB spectra or correlations
pars.Want_CMB = False 
#NonLinear_both or NonLinear_lens will use non-linear corrections
pars.NonLinear = model.NonLinear_both
#Set up W(z) window functions, later labelled W1, W2. Gaussian here.
pars.SourceWindows = [
    GaussianSourceWindow(redshift=0.17, source_type='counts', bias=1.2, sigma=0.04, dlog10Ndm=-0.2),
    GaussianSourceWindow(redshift=0.5, source_type='lensing', sigma=0.07)]

results = camb.get_results(pars)
cls = results.get_source_cls_dict()

#Note that P is CMB lensing, as a deflection angle power (i.e. PxP is [l(l+1)]^2C_l\phi\phi/2\pi)
#lensing window functions are for kappa (and counts for the fractional angular number density)
ls=  np.arange(2, lmax+1)
for spectrum in ['W1xW1','W2xW2','W1xW2',"PxW1", "PxW2"]:
    plt.loglog(ls, cls[spectrum][2:lmax+1], label=spectrum)
plt.xlabel(r'$\ell$')
plt.ylabel(r'$\ell(\ell+1)C_\ell/2\pi$')
plt.legend();

#You can also use window functions from numerical table of W(z). It must be well enough sampled to interpolate nicely.
#e.g. reproduce Gaussian this way..
zs = np.arange(0, 0.5, 0.02)
W = np.exp(-(zs - 0.17) ** 2 / 2 / 0.04 ** 2) / np.sqrt(2 * np.pi) / 0.04
pars.SourceWindows = [SplinedSourceWindow(bias=1.2, dlog10Ndm=-0.2, z=zs, W=W)]
results = camb.get_results(pars)
cls2 = results.get_cmb_unlensed_scalar_array_dict()
plt.loglog(ls, cls["W1xW1"][2:lmax+1], label=spectrum)
plt.loglog(ls, cls2["W1xW1"][2:lmax+1],ls='--')
plt.xlabel(r'$\ell$')
plt.ylabel(r'$\ell(\ell+1)C_\ell/2\pi$');

#Sources can include various terms using these options (line_xx refers to 21cm)
print(pars.SourceTerms)

class: <SourceTermParams>
 limber_windows = True
 limber_phi_lmin = 100
 counts_density = True
 counts_redshift = True
 counts_lensing = False
 counts_velocity = True
 counts_radial = False
 counts_timedelay = True
 counts_ISW = True
 counts_potential = True
 counts_evolve = False
 line_phot_dipole = False
 line_phot_quadrupole = False
 line_basic = True
 line_distortions = True
 line_extra = False
 line_reionization = False
 use_21cm_mK = True
 

#Results above include redshift distortions but not magnification bias (counts_lensing). 
#Try turning off redshift distortions:
pars.SourceTerms.counts_redshift = False
results = camb.get_results(pars)
cls3 = results.get_source_cls_dict()

plt.loglog(ls, cls["W1xW1"][2:lmax+1])
plt.loglog(ls, cls3["W1xW1"][2:lmax+1])
plt.legend(['With redshift distortions', 'Without'])
plt.xlabel(r'$\ell$')
plt.ylabel(r'$\ell(\ell+1)C_\ell/2\pi$')
plt.xlim(2,lmax);

# For number counts you can give a redshift-dependent bias (the underlying code supports general b(z,k))
# toy model for single-bin LSST/Vera Rubin [using numbers from 1705.02332]
z0=0.311
zs = np.arange(0, 10, 0.02)
W = np.exp(-zs/z0)*(zs/z0)**2/2/z0
bias = 1 + 0.84*zs
pars.SourceWindows = [SplinedSourceWindow(dlog10Ndm=0, z=zs, W=W, bias_z = bias)]
lmax=3000
pars.set_for_lmax(lmax, lens_potential_accuracy=5)
results = camb.get_results(pars)

#e.g. plot the cross-correlation with CMB lensing
cls = results.get_cmb_unlensed_scalar_array_dict()

nbar = 40/(np.pi/180/60)**2 # Poission noise
ls= np.arange(2,lmax+1)
Dnoise = 1/nbar*ls*(ls+1)/2/np.pi

correlation=cls["PxW1"][2:lmax+1]/np.sqrt(cls["PxP"][2:lmax+1]*(cls["W1xW1"][2:lmax+1]+Dnoise))
plt.plot(np.arange(2,lmax+1), correlation)
plt.xlabel(r'$L$')
plt.ylabel('correlation')
plt.xlim(2,lmax)
plt.title('CMB lensing - LSST correlation (single redshift bin)');

#Let's look at some non-standard primordial power spectrum, e.g. with wavepacket oscillation

#Define our custom  power spectrum function (here power law with one wavepacket)
def PK(k, As, ns, amp, freq, wid, centre, phase):
    return As*(k/0.05)**(ns-1)*(1+ np.sin(phase+k*freq)*amp*np.exp(-(k-centre)**2/wid**2))

#Check how this looks compared to power law
ks = np.linspace(0.02,1,1000)
pk1 = 2e-9*(ks/0.05)**(0.96-1)
pk2 = PK(ks,2e-9, 0.96,0.0599, 280, 0.08, 0.2,0)
plt.semilogx(ks,pk1)
plt.semilogx(ks,pk2)
plt.ylabel('$P(k)$')
plt.xlabel(r'$k\, {\rm Mpc}$')
plt.legend(['Power law','Custom'])
plt.title('Scalar initial power spectrum');

#Now compute C_l and compare
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, tau=0.06)
lmax=2500
pars.set_for_lmax(lmax,lens_potential_accuracy=1)

#For comparison, standard power law
pars.InitPower.set_params(As=2e-9, ns=0.96)
results = camb.get_results(pars)
cl_unlensed=results.get_unlensed_scalar_cls(CMB_unit ='muK')
cl=results.get_lensed_scalar_cls(CMB_unit ='muK')

#Now get custom spectrum (effective_ns_for_nonlinear is used for halofit if required)
pars.set_initial_power_function(PK, args=(2e-9, 0.96,0.0599, 280, 0.08, 0.2,0), 
                                effective_ns_for_nonlinear=0.96)

results2 = camb.get_results(pars)
cl2=results2.get_lensed_scalar_cls(CMB_unit ='muK')

ls = np.arange(2,lmax)
plt.plot(ls,(cl2[2:lmax,0]-cl[2:lmax,0]))
plt.xlabel(r'$\ell$');
plt.ylabel(r'$\ell(\ell+1)\Delta C_\ell/2\pi\, [\mu K^2]$')
plt.title(r'$C_\ell$ difference to power law');

#Note that if you have sharp features or fine oscillations, you may need 
#increase accuracy to sample them well. e.g. let's try increasing the frequency

#Default accuracy
pars.Accuracy.lSampleBoost = 1
pars.Accuracy.IntkAccuracyBoost =1
pars.Accuracy.SourcekAccuracyBoost =1

freq = 1000
ks = np.linspace(0.02,1,1000)
pk1 = 2e-9*(ks/0.05)**(0.96-1)
pk2 = PK(ks,2e-9, 0.96,0.0599,freq, 0.08, 0.2,0)
plt.semilogx(ks,pk1)
plt.semilogx(ks,pk2)
plt.ylabel('$P(k)$')
plt.xlabel(r'$k\, {\rm Mpc}$');
plt.title('Scalar power spectrum')
plt.figure()

pars.set_initial_power_function(PK, args=(2e-9, 0.96,0.0599, freq, 0.08, 0.2,0),
                               effective_ns_for_nonlinear=0.96)

results2 = camb.get_results(pars)
cl_unlensed2=results2.get_unlensed_scalar_cls(CMB_unit ='muK')

#need to increase default sampling in ell to see features smaller than peaks reliably
pars.Accuracy.lSampleBoost = 2
#may also need to sample k more densely when computing C_l from P(k)
pars.Accuracy.IntkAccuracyBoost = 2 

results3 = camb.get_results(pars)
cl_unlensed3=results3.get_unlensed_scalar_cls(CMB_unit ='muK')
cl3=results3.get_lensed_scalar_cls(CMB_unit ='muK')

ls = np.arange(2,lmax)
plt.plot(ls,(cl_unlensed2[2:lmax,0]/cl_unlensed[2:lmax,0]-1))
plt.plot(ls,(cl_unlensed3[2:lmax,0]/cl_unlensed[2:lmax,0]-1))

plt.xlabel(r'$\ell$')
plt.ylabel(r'$\Delta C_\ell/C_\ell$')
plt.title(r'Fractional $C_\ell$ difference to power law')
plt.legend(['Default accuracy','Boosted accuracy']);

#Note that lensing washes out small features on small scales
plt.plot(ls,(cl_unlensed3[2:lmax,0]/cl_unlensed[2:lmax,0]-1))
plt.plot(ls,(cl3[2:lmax,0]/cl[2:lmax,0]-1))
plt.xlabel(r'$\ell$')
plt.ylabel(r'$\Delta C_\ell/C_\ell$')
plt.legend(['Unlensed','Lensed'])
plt.title(r'Fractional $C_\ell$ difference to power law');

#Now look at the (small!) effect of neutrino mass splittings on the matter power spectrum (in linear theory)
#The "neutrino_hierarchy" parameter uses a two eigenstate approximation to the full hierarchy (which is very accurate for cosmology)
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.11, neutrino_hierarchy='normal')
pars.InitPower.set_params(ns=0.965)
pars.set_matter_power(redshifts=[0.], kmax=2.0)
results = camb.get_results(pars)
kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=2, npoints = 200)

pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.11, neutrino_hierarchy='inverted')
results = camb.get_results(pars)
kh2, z2, pk2 = results.get_matter_power_spectrum(minkh=1e-4, maxkh=2, npoints = 200)

plt.semilogx(kh, pk2[0,:]/pk[0,:]-1)
plt.ylabel('$\Delta$ PK/PK')
plt.xlabel(r'$k\, [h \,\rm{Mpc}^{-1}]$')
plt.title(r'Normal vs Inverted for $\sum m_\nu=0.11 \rm{eV}$');

#Matter power functions can get other variables than the total density.
#For example look at the relative baryon-CDM velocity by using the power spectrum of
#model.Transfer_vel_baryon_cdm

kmax=10 
k_per_logint = 30
zs = [200, 500,800, 1090]

pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(ns=0.965)
pars.WantTransfer = True
pars.set_matter_power(redshifts=zs, kmax=kmax, k_per_logint=k_per_logint, silent=True)
results = camb.get_results(pars)


PKint = results.get_matter_power_interpolator(nonlinear=False, 
    hubble_units=False, k_hunit=False, 
    var1=model.Transfer_vel_baryon_cdm, var2=model.Transfer_vel_baryon_cdm)

#Make plot like Fig 1 of https://arxiv.org/abs/1005.2416
#showing contribution to the relative velocity variance per log k

plt.figure(figsize=(8,5))
ks = np.logspace(-3, np.log10(kmax), 300)
for i, z in enumerate(zs):
     curlyP = PKint.P(z,ks)*ks**3/(2*np.pi**2)
     plt.semilogx(ks, curlyP)
     rms = np.sqrt(curlyP[1:-1].dot((ks[2:]-ks[:-2])/2/ks[1:-1]))
     print('rms velocity at z=%s: %.3g'%(z,rms), '(%.3g km/s)'%(rms*299792))
    
plt.xlim([1e-3,kmax])
plt.xlabel('k Mpc', fontsize=16)
plt.ylabel(r'$\mathcal{P}_v(k)$', fontsize=16)
plt.legend(['$z = %s$'%z for z in zs]);

rms velocity at z=200: 2.04e-05 (6.11 km/s)
rms velocity at z=500: 5.03e-05 (15.1 km/s)
rms velocity at z=800: 7.89e-05 (23.6 km/s)
rms velocity at z=1090: 9.67e-05 (29 km/s)

#You can also get the matter transfer functions
#These are synchronous gauge and normalized to unit primordial curvature perturbation
#The values stored in the array are quantities like Delta_x/k^2, and hence
#are nearly independent of k on large scales. 
#Indices in the transfer_data array are the variable type, the k index, and the redshift index

pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(ns=0.965)
pars.set_matter_power(redshifts=[0], kmax=kmax)
results= camb.get_results(pars)

trans = results.get_matter_transfer_data()
#get kh - the values of k/h at which they are calculated
kh = trans.transfer_data[0,:,0]
#transfer functions for different variables, e.g. CDM density and the Weyl potential
#CDM perturbations have grown, Weyl is O(1) of primordial value on large scales
delta = trans.transfer_data[model.Transfer_cdm-1,:,0]
W = trans.transfer_data[model.Transfer_Weyl-1,:,0]
plt.plot(kh,delta)
plt.loglog(kh,-W)
plt.xlabel(r'$k/h\, [\rm Mpc]^{-1}$', fontsize=16);
plt.title('Matter transfer functions')
plt.legend([r'$\Delta_c/k^2$',r'Weyl potential $\Psi$'], fontsize=14);

#Check we can get the matter power spectrum from the transfer function as expected
k = kh*results.Params.h
transfer = trans.transfer_data[model.Transfer_tot-1,:,0]
primordial_PK = results.Params.scalar_power(k)
matter_power = primordial_PK*transfer**2*k**4 / (k**3/(2*np.pi**2))

#compare with CAMB's explicit output for the matter power spectrum
kh2,zs,PK = results.get_linear_matter_power_spectrum(hubble_units=False)

plt.loglog(kh,matter_power)
plt.loglog(kh, PK[0,:]);
plt.xlabel(r'$k\, [h Mpc^{-1}]$');

# It is also possible to get the real-space linear perturbation variance in spheres, i.e. \sigma_R.
# Using R=8 and hubble_units=T would return the standard definition of \sigma_8
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(ns=0.965)
pars.set_matter_power(redshifts=[0,1], kmax=kmax)
results= camb.get_results(pars)
R, z, sigma_R = results.get_sigmaR(R=np.arange(1,20,0.5), 
                                 hubble_units=False, return_R_z=True)
plt.plot(R, sigma_R[1,:], label='z = %s'%z[1])
plt.plot(R, sigma_R[0,:], label='z = %s'%z[0])
plt.ylabel(r'$\sigma_R$', fontsize=14)
plt.xlabel(r'$R/{\rm Mpc}$', fontsize=14)
plt.legend()
# To get the non-linear scale where sigma_R=1, e.g. at z=0
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.optimize import brentq
sR = InterpolatedUnivariateSpline(R, sigma_R[-1,:]-1)
R_nonlin = brentq(sR, R[0], R[-1])
print(r'R giving \sigma_R=1 at z=0 is at R=%.2f Mpc (or %.2f Mpc/h)'
                         %(R_nonlin, R_nonlin * pars.h))

Note: redshifts have been re-sorted (earliest first)
R giving \sigma_R=1 at z=0 is at R=8.51 Mpc (or 5.74 Mpc/h)

# the above is for the default density variable delta_tot, without including
# massive neutrinos the result would be very slightly different
sigma_R_nonu = results.get_sigmaR(R=np.arange(1,20,0.5), 
                                 var1='delta_nonu', var2='delta_nonu',
                                 hubble_units=False)
plt.plot(R, sigma_R_nonu[0,:]/ sigma_R[0,:]-1, label='z = %s'%z[0])
plt.plot(R, sigma_R_nonu[1,:]/ sigma_R[1,:]-1, label='z = %s'%z[1])
plt.ylabel(r'$\Delta\sigma_R/\sigma_R$', fontsize=14)
plt.xlabel(r'$R/{\rm Mpc}$', fontsize=14)
plt.legend();

# results for power spectra using different initial power spectra 
# can be computed without  recomputing the transfer functions

pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, redshifts=[0], kmax=5, 
                       As=2e-9, ns=0.96, halofit_version='mead')
results= camb.get_transfer_functions(pars)
kref,_,PKref = results.get_linear_matter_power_spectrum(hubble_units=False, k_hunit=False)
kref,_,PKnl = results.get_nonlinear_matter_power_spectrum(hubble_units=False, k_hunit=False)

ns_arr = np.arange(0.94, 0.981, 0.01)
fig, axs = plt.subplots(1,2, figsize=(14,5), sharey=True)

for ns in ns_arr:
    results.Params.InitPower.set_params(As=2e-9, ns=ns)
    results.calc_power_spectra()
    k,_,PK = results.get_linear_matter_power_spectrum(hubble_units=False, k_hunit=False)
    np.testing.assert_allclose(k,kref)
    axs[0].semilogx(k, PK[0,:]/PKref[0,:]-1);
    k,_,PK = results.get_nonlinear_matter_power_spectrum(hubble_units=False, k_hunit=False)
    axs[1].semilogx(k, PK[0,:]/PKnl[0,:]-1);

for ax, t in zip(axs,['Linear', 'Nonlinear']):
    ax.set_title('%s spectrum'%t)
    ax.set_xlim(1e-3,5)
    #ax.set_ylim([-0.1,0.1])
    ax.set_xlabel(r'$k\, [Mpc^{-1}]$');
plt.legend(['$n_s = %.2f$'%ns for ns in ns_arr], ncol=2, loc='lower left');
axs[0].set_ylabel(r'$\Delta P(k)/P(k)$', fontsize=16);

# The non-linear model parameters can also be varied without recomputing the transfer functions
# eg. look at the effect of the HMcode baryonic feedback parameter
results= camb.get_transfer_functions(pars)

kref,_,PKnl = results.get_nonlinear_matter_power_spectrum(hubble_units=False, k_hunit=False)
feedbacks= np.arange(2, 4.1, 0.5)
for baryfeed in feedbacks:
    results.Params.NonLinearModel.set_params(halofit_version='mead', HMCode_A_baryon=baryfeed)
    results.calc_power_spectra()
    k,_,PK = results.get_nonlinear_matter_power_spectrum(hubble_units=False, k_hunit=False)
    np.testing.assert_allclose(k,kref)
    plt.semilogx(k, PK[0,:]/PKnl[0,:]-1)
    plt.xlim(1e-2,5)
plt.ylabel(r'$\Delta P(k)/P(k)$', fontsize=16)
plt.legend([r'$A_{\rm baryon} = %.2f$'%b for b in feedbacks]);

# However non-linear lensing or other sources, or non-linearly lensed CMB requires 
# recalculation from the time transfer functions. This is a bit slower but faster than recomputing everything.
# e.g. look at the impact of the baryon feedback parameter on the lensing potential power spectrum

pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, ns=0.965, lens_potential_accuracy=1, lmax=3000)
results= camb.get_transfer_functions(pars, only_time_sources=True)
results.calc_power_spectra()
ref = results.get_lens_potential_cls()[:,0]

feedbacks= np.arange(2, 4.1, 0.5)
for baryfeed in feedbacks:
    results.Params.NonLinearModel.set_params(halofit_version='mead', HMCode_A_baryon=baryfeed)
    results.calc_power_spectra()
    CL = results.get_lens_potential_cls()[:,0]
    plt.semilogx(np.arange(2,CL.shape[0]),CL[2:]/ref[2:]-1);

plt.legend([r'$A_{\rm baryon} = %.2f$'%b for b in feedbacks])
plt.ylabel(r'$\Delta C_L/C_L$', fontsize=16)
plt.xlabel('$L$')
plt.title('Impact of baryonic feedback on the lensing potential');

# CAMB has a basic scalar field quintessence dark energy model
# EarlyQuintessence is a specific example implementation that implements early dark energy
# (axion-like, as arXiv:1908.06995) with potential V(\phi) = m^2f^2 (1 - cos(\phi/f))^2 + \Lambda
# AxionEffectiveFluid is an approximate model that does not evolve the quintessence equations
# To implement other models you'd need to make your own new subclass.
# Note that this is not as well tested as most of the rest of CAMB.

# Use n=3 and keep theta_* angular distance parameter fixed to roughly fit CMB data
thetastar= 0.01044341764253 
n=3

fig, axs = plt.subplots(3,1, figsize=(10,8))
zs = np.logspace(1,5,500)
pars = camb.set_params( ombh2=0.022, omch2=0.122, thetastar=thetastar) #
results = camb.get_results(pars)
print('LCDM: h0=', results.Params.H0)
cl_LCDM = results.get_lensed_scalar_cls(CMB_unit='muK')
axs[1].plot(cl_LCDM[2:,0])

# Set dark energy fraction 0.1 at z=1e4
pars = camb.set_params( ombh2=0.022, omch2=0.122, 
                       w_n=(n-1.)/(n+1.),  theta_i=np.pi/2, zc = 1e4, fde_zc = 0.1,
                       dark_energy_model='AxionEffectiveFluid', thetastar=thetastar) #
results = camb.get_results(pars)
print('AxionEffectiveFluid: h0 = ', results.Params.H0)
axs[0].semilogx(zs,results.get_Omega('de',z=zs))
cl = results.get_lensed_scalar_cls(CMB_unit='muK')
axs[0].set_ylabel(r'$\Omega_{\rm de}$')
axs[2].plot(cl[2:,0]/cl_LCDM[2:,0]-1)

pars = camb.set_params( ombh2=0.022, omch2=0.122, 
                       m=8e-53, f =0.05,n=n,  theta_i=3.1,use_zc = True, zc = 1e4, fde_zc = 0.1,
                       dark_energy_model='EarlyQuintessence', thetastar=thetastar) #
results = camb.get_results(pars)
print('EarlyQuintessence: h0 = ', results.Params.H0)
cl = results.get_lensed_scalar_cls(CMB_unit='muK')
axs[0].semilogx(zs,results.get_Omega('de',z=zs))
axs[0].set_xlabel(r'$z$')
axs[1].plot(cl[2:,0])
axs[2].plot(cl[2:,0]/cl_LCDM[2:,0]-1)
axs[1].set_ylabel(r'$D_l [\mu {\rm K}^2]$')
axs[2].set_ylabel(r'$\Delta D_l [\mu {\rm K}^2]$')
axs[1].set_xlabel(r'$\ell$')
axs[2].set_xlabel(r'$\ell$')
plt.tight_layout()

print('m = ',results.Params.DarkEnergy.m,'f = ',results.Params.DarkEnergy.f )

LCDM: h0= 67.40599027170882
AxionEffectiveFluid: h0 =  78.296070152354
EarlyQuintessence: h0 =  73.21865482559397
m =  2.1483075882436665e-54 f =  0.07349231434987724

# You can also calculate CMB correlation functions
# (the correlations module also has useful functions for the inverse transform,
# CMB lensing and derivatives of the lensed spectra - see docs)
from camb import correlations

pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, As=2e-9, ns=0.96)
pars.set_for_lmax(4000, lens_potential_accuracy=1)
results = camb.get_results(pars)
lensed_cl = results.get_lensed_scalar_cls()
corrs, xvals, weights = correlations.gauss_legendre_correlation(lensed_cl)

r=np.arccos(xvals)*180/np.pi # sampled theta values in degrees
fig, axs= plt.subplots(2,2, figsize=(12,8))
for ix, (cap, ax) in enumerate(zip(['TT','+','-',r'\times'], axs.reshape(4))):
    ax.plot(r, corrs[:,ix])
    ax.axhline(0,color='k')
    ax.set_xlim([0,10])
    ax.set_xlabel(r'$\theta$')
    ax.set_ylabel(r'$\zeta_{%s}(\theta)$'%(cap), fontsize=14)
    ax.set_xlim(0,4)
plt.suptitle('Correlation functions');        

# For partially-delensed or Alens-scaled spectra, power spectra can be computed using 
# a custom or scaled lensing potential power spectrum. There's a pure-python 
# interface in the correlation module, or can use the result object functions 
# (faster).
# Here just plot results for scaled lensing spectrum, can use 
# get_lensed_cls_with_spectrum to calculate lensed spectrum with specific 
# lensing power if needed. For BB the scaling is fairly linear, but less so for
# other spectra.

pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, As=2e-9, ns=0.96)
pars.set_for_lmax(3500, lens_potential_accuracy=1)
results = camb.get_results(pars)
for Alens in [1,0.9, 0.5, 0.1]:
    plt.plot(np.arange(2,2501),
             results.get_partially_lensed_cls(Alens, lmax=2500)[2:,2],
            label='$A_L = %s$'%Alens)
plt.ylabel(r'$C_l^{BB}$')
plt.xlabel('$l$')   
plt.xlim([2,2500])
plt.legend();

# For CMB lensing reconstruction, the non-perturbative response functions depend on the 
# `gradient' C_l, which can be calculated in the flat sky approximation using CAMB
# On large scales similar to lensed CMB spectrum, but important difference esp for TT on small scales
# see arXiv:1906.08760, 1101.2234 
c_lensed = results.get_lensed_scalar_cls()
c_lens_response = results.get_lensed_gradient_cls()
plt.plot(np.arange(2,3501), c_lens_response[2:3501,0]/c_lensed[2:3501,0]-1)
plt.ylabel(r'$\Delta C_l/C_l$')
plt.xlabel('$l$')  
plt.xlim(2,3500)
plt.title(r'Lensing gradient response vs lensed $C_l$');







CAMB also supports dark-age 21cm (see astro-ph/0702600)

 







 





# Get 21cm transfer functions, similar to Fig 4 of astro-ph/0702600
pars=camb.set_params(cosmomc_theta=0.0104, ombh2= 0.022, omch2= 0.12, tau = 0.053,
                     As= 2.08e-09, ns= 0.968, num_massive_neutrinos=1)

pars.Do21cm  =True

pars.Evolve_delta_xe =True # ionization fraction perturbations
pars.Evolve_baryon_cs  = True # accurate baryon perturbations
pars.WantCls =False
pars.SourceTerms.use_21cm_mK = False # Use dimensionless rather than mK units
redshifts=[50]
pars.set_matter_power(kmax=1000, redshifts=redshifts)

# get transfer functions
results= camb.get_results(pars)
trans = results.get_matter_transfer_data()
#get kh - the values of k/h at which they are calculated
k = trans.transfer_data[0,:,0]* results.Params.h

for zix, z in enumerate(redshifts,):
    # see https://camb.readthedocs.io/en/latest/transfer_variables.html

    delta_c = trans.transfer_data[model.Transfer_cdm-1,:,zix]
    delta_b = trans.transfer_data[model.Transfer_b-1,:,zix]
    T_g = trans.transfer_data[model.Transfer_Tmat-1,:,zix]
    mono = trans.transfer_data[model.Transfer_monopole-1,:,zix]

    plt.loglog(k,mono*k**2, color='k', ls='-', lw=2)
    plt.plot(k,delta_b*k**2, color='gray', ls='-')
    plt.plot(k,T_g*k**2, color='b', ls='--')
    plt.plot(k,delta_c*k**2, color='C0', ls='-.')

plt.xlabel(r'$k\, \rm Mpc$', fontsize=16);
plt.title(f'Transfer functions, z={redshifts}')
plt.legend([r'21cm monopole', r'$\Delta_b$',r'$\Delta_{T_g}$', r'$\Delta_c$'], fontsize=14)
plt.ylabel(r'$T_\chi$')
plt.xlim(1e-4, 1e3)
plt.ylim(1e-3, 1e4);


 







 

  

 

       

 

 

 

 



CAMB also supports dark-age 21cm (see astro-ph/0702600)

# Get 21cm transfer functions, similar to Fig 4 of astro-ph/0702600
pars=camb.set_params(cosmomc_theta=0.0104, ombh2= 0.022, omch2= 0.12, tau = 0.053,
                     As= 2.08e-09, ns= 0.968, num_massive_neutrinos=1)

pars.Do21cm  =True

pars.Evolve_delta_xe =True # ionization fraction perturbations
pars.Evolve_baryon_cs  = True # accurate baryon perturbations
pars.WantCls =False
pars.SourceTerms.use_21cm_mK = False # Use dimensionless rather than mK units
redshifts=[50]
pars.set_matter_power(kmax=1000, redshifts=redshifts)

# get transfer functions
results= camb.get_results(pars)
trans = results.get_matter_transfer_data()
#get kh - the values of k/h at which they are calculated
k = trans.transfer_data[0,:,0]* results.Params.h

for zix, z in enumerate(redshifts,):
    # see https://camb.readthedocs.io/en/latest/transfer_variables.html

    delta_c = trans.transfer_data[model.Transfer_cdm-1,:,zix]
    delta_b = trans.transfer_data[model.Transfer_b-1,:,zix]
    T_g = trans.transfer_data[model.Transfer_Tmat-1,:,zix]
    mono = trans.transfer_data[model.Transfer_monopole-1,:,zix]

    plt.loglog(k,mono*k**2, color='k', ls='-', lw=2)   
    plt.plot(k,delta_b*k**2, color='gray', ls='-')
    plt.plot(k,T_g*k**2, color='b', ls='--')
    plt.plot(k,delta_c*k**2, color='C0', ls='-.')
    
plt.xlabel(r'$k\, \rm Mpc$', fontsize=16);
plt.title(f'Transfer functions, z={redshifts}')
plt.legend([r'21cm monopole', r'$\Delta_b$',r'$\Delta_{T_g}$', r'$\Delta_c$'], fontsize=14)
plt.ylabel(r'$T_\chi$')
plt.xlim(1e-4, 1e3)
plt.ylim(1e-3, 1e4);

Animations

Going back to standard CMB, here is an example of how to to make an animation of the evolution of the transfer functions.

# A potential issue here is that with large dynamic range, you may wish to plot modes where
# k*eta << 1 (long way outside horizon). Evolution is not normally calculated in the code a long
# way outside the horizon, starting instead with the series solution. So for results which are numerically
# unstable, you may need to replace the numerical result with the series result.

# Use widget to see animation in notebook rather than just saving (with "pip install ipympl")
# %matplotlib widget 

pars = camb.read_ini('https://raw.githubusercontent.com/cmbant/CAMB/master/inifiles/planck_2018.ini')
data= camb.get_results(pars)

# get ranges to plot, evolving until recombination
zstar=data.get_derived_params()['zstar']
etastar=data.conformal_time(zstar)
# stop time evolution (eta) at recombination
eta = np.logspace(-1.7,np.log10(etastar),400)
# wide range of k
k=10**np.linspace(-3,1,1200)

# get some background quantities
z = data.redshift_at_conformal_time(eta)
back_ev = data.get_background_time_evolution(eta, ['x_e', 'visibility'])
x_e=back_ev['x_e']
vis=back_ev['visibility']
adotoa = data.h_of_z(z)/(1+z)  # comoving Hubble
# ratio of matter to radiation (neglecting neutrino mass)
rho=data.get_background_densities(1/(1+z))
Rmat=(rho['baryon']+rho['cdm'])/(rho['photon']+rho['neutrino']+rho['nu'])

# some quantities needed for superhorizon series solution -see series solutions in notes at
# https://cosmologist.info/notes/CAMB.pdf
grhonu=data.grhormass[0]+data.grhornomass
Rv=grhonu/(grhonu+data.grhog)
omega = (data.grhob+data.grhoc)/np.sqrt(3*(data.grhog+grhonu))
# initial value of super-horizon Psi for unit curvature
Psi_init = -10/(4*Rv+15)

# make evolution plot in Newtonian gauge. Use symbolic package to get right gauge-invariant quantities.
import camb.symbolic as cs

ev=data.get_time_evolution(k, eta, [cs.Psi_N, # Newtonian potential \Psi
                                    'delta_photon', # sychronous gauge photon perturbation
                                    cs.make_frame_invariant(cs.Delta_g, 'Newtonian'), #photon perturbation
                                    cs.make_frame_invariant(cs.v_b, 'Newtonian') # baryon velocity
                                   ], lAccuracyBoost=4)
Psi = ev[:,:,0]
Delta_g = ev[:,:,1]
Delta_g_N=ev[:,:,2]
v_b_N = ev[:,:,3]

# Now replace Psi results well outside horizon where numerically unstable with series solution
# Note Newtonian-gauge Delta_g does not start at zero, so is also unstable
for i, etai in enumerate(eta):
    for j, kj in enumerate(k):
        if etai*kj < 0.1 and etai*omega < 0.1:
            Psi[j,i]=Psi_init - 25/8*omega*etai *(8*Rv-3)/(4*Rv+15)/(2*Rv+15)
            # use Delta_g_N = Delta_g - 4H sigma/k (with adiabatic series result for sigma)
            Delta_g_N[j,i] = Delta_g[j,i] - 4*adotoa[i]*( Psi_init/2*etai - 15*omega*etai**2/8*(4*Rv-5)/(4*Rv+15)/(2*Rv+15))

import math
from matplotlib import animation

def latex_sci_format(num, dec=3):
    return "{:.{dec}f}\\times 10^{{{}}}".format(num / (10 ** int(math.floor(math.log10(abs(num))))), int(math.floor(math.log10(abs(num)))), dec=dec)

fig, ax = plt.subplots()

plot_vars = {r'$\Psi$':Psi, 
             r'$\frac{1}{4}\Delta^{\rm Newt}_\gamma +\Psi$': Delta_g_N/4+Psi, 
             r'$\frac{1}{\sqrt{3}}v_b^{\rm Newt}$': v_b_N/np.sqrt(3) }

anim_fps=8

lines=[]
for lab in plot_vars.keys():
    line, = ax.semilogx([], [], label=lab)
    lines.append(line)
    
ax.axhline(0,color='k')
plt.legend(loc='upper left')

ax.set_xlim(k[0], k[-1])
ax.set_ylim(-1, 1)
ax.set_xlabel(r'$k\, {\rm Mpc}$')

# Define the update function
def update(i):
    ax.set_title( r'$z= %s, \eta/{\rm Mpc}= %.2f, x_e = %.2f, \rho_m/\rho_{\rm rad}=%.2f$' 
                  % (latex_sci_format(z[i]), eta[i], x_e[i], Rmat[i]))
    for line, var in zip(lines,plot_vars.values()):
        line.set_data(k, var[:, i])
    return lines

# Create the animation
transfer_anim =animation.FuncAnimation(fig, update, frames=range(len(eta)),  blit=True)
writer = animation.writers['ffmpeg'](fps=anim_fps, bitrate=-1)

# can save like this
transfer_anim.save(r'Newtonian_transfer_evolve_monopole_phi_vb.mp4', writer=writer, dpi=240)

See the output mp4 file here or play below:

# or play inline using:
# from IPython.display import display, Video
# video = transfer_anim.to_html5_video()
# display(Video(data=video, embed=True, width = 600))

# Animate the matter transfer and potential evolution up to today
# now use log scale and synchronous gauge as more relevant for the late-time matter

eta=np.logspace(-1.7,np.log10(data.tau0),500)[:-1]
k=10**np.linspace(-3,2,2000)
z = data.redshift_at_conformal_time(eta)
import camb.symbolic as cs
ev=data.get_time_evolution(k, eta, [cs.Psi_N,
                                    'delta_cdm', # sychronous gauge cdm perturbation
                                    'delta_photon', 
                                    'delta_baryon'], lAccuracyBoost=4)

Psi = ev[:,:,0]
Delta_c = ev[:,:,1]
Delta_g=ev[:,:,2]
Delta_b=ev[:,:,3]

x_e = data.get_background_time_evolution(eta, ['x_e'])['x_e']
adotoa = data.h_of_z(z)/(1+z)  # comoving Hubble
rho=data.get_background_densities(1/(1+z))
Rmat=(rho['baryon']+rho['cdm'])/(rho['photon']+rho['neutrino']+rho['nu'])
zstar=data.get_derived_params()['zstar']

# series for super-horizon Psi
grhonu=data.grhormass[0]+data.grhornomass
Rv=grhonu/(grhonu+data.grhog)
omega = (data.grhob+data.grhoc)/np.sqrt(3*(data.grhog+grhonu))
# initial value of super-horizon Psi for unit curvature and no isocurvature
Psi_init = -10/(4*Rv+15)
for i, etai in enumerate(eta):
    for j, kj in enumerate(k):
        if etai*kj < 0.1 and etai*omega < 0.1:
            Psi[j,i]=Psi_init - 25/8*omega*etai *(8*Rv-3)/(4*Rv+15)/(2*Rv+15)

fig, ax = plt.subplots()
plot_vars = {r'$\Delta_c$': Delta_c, r'$\Delta_b$': Delta_b,
             r'$\Delta_g$': Delta_g,r'$\Psi$': Psi}             
anim_fps=8
lines=[]
for lab in plot_vars.keys():
    line, = ax.loglog([], [], label=lab)
    lines.append(line)

plt.axhline(0, ls='-', color='k')
plt.legend(loc='upper left')
ax.set_xlim(k[0], k[-1])
ax.set_xlabel(r'$k\, {\rm Mpc}$')

ax.set_ylim(1e-4, 2e5)
lines.append(ax.loglog([], [],ls='--', color='C1')[0])
lines.append(ax.loglog([], [],ls='--', color='C2')[0])
lines.append(ax.loglog([], [],ls='--', color='C3')[0])
    
# Define the update function
def update(i):
    ax.set_title( r'$z= %s, \eta/{\rm Mpc}= %.2f, x_e = %.2f, \rho_m/\rho_{\rm rad}=%.2f$' 
                  % (latex_sci_format(z[i],2), eta[i], x_e[i], Rmat[i]), fontsize=10)

    for line, var in zip(lines,plot_vars.values()):
        line.set_data(k, var[:, i])
    lines[-3].set_data(k,-Delta_b[:, i])
    lines[-2].set_data(k,-Delta_g[:, i])
    lines[-1].set_data(k,-Psi[:, i])
        
    if z[i] < zstar: # hide photons after recombination as irrelevant for matter
        lines[2].set_data([],[]) 
        lines[-2].set_data([],[])      
    return lines

# Create the animation
transfer_anim =animation.FuncAnimation(fig, update, frames=range(len(eta)),  blit=True, repeat = False)
writer = animation.writers['ffmpeg'](fps=anim_fps, bitrate=-1)
transfer_anim.save(r'matter_transfer_evolve_sync.mp4', writer=writer, dpi=240)

C:\Users\Antony\AppData\Local\Temp\ipykernel_87596\455262666.py:12: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
  ax.set_xlim(k[0], k[-1])

This the video produced of the synchronous transfer function evolution (dashed lines are negative):